rm(list=ls())
gc()
setwd("/hpc/group/pbenfeylab/CheWei/")
| used | (Mb) | gc trigger | (Mb) | max used | (Mb) | |
|---|---|---|---|---|---|---|
| Ncells | 625275 | 33.4 | 1361470 | 72.8 | 1209572 | 64.6 |
| Vcells | 1159072 | 8.9 | 8388608 | 64.0 | 1802279 | 13.8 |
as.numeric(system("awk '/MemFree/ {print $2}' /proc/meminfo", intern=TRUE))
suppressMessages(library(Seurat))
suppressMessages(library(tricycle))
suppressMessages(library(cowplot))
suppressMessages(library(ggplot2))
suppressMessages(library(scattermore))
suppressMessages(library(scater))
suppressMessages(library(cowplot))
suppressMessages(library(RColorBrewer))
suppressMessages(library(grid))
suppressMessages(library(gplots))
suppressMessages(library(tidyverse))
prep <- function(seu){
seu$time.anno.Li.crude <- gsub("_.*$","",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Columella","Distal Columella_Columella",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Lateral Root Cap","Distal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Columella","Proximal Columella_Columella",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Lateral Root Cap","Proximal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("^.*_","",seu$celltype.anno.Li.crude)
return(seu)
}
plot_anno <- function(rc.integrated){
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li.crude <- factor(rc.integrated$celltype.anno.Li.crude, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]
p1 <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li.crude", cols=color)
p2 <- DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno.Li.crude", order = c("Distal Columella","Proximal Columella","Distal Lateral Root Cap","Proximal Lateral Root Cap","Maturation","Elongation", "Transition Domain", "Proliferation Domain"),
cols = c('#F7E7B0','#FFC400','#2B871F','#005E3B',"#deebf7", "#3182bd", '#fee0d2','#de2d26'))
p3 <- DimPlot(rc.integrated, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral"))
p4 <- DimPlot(rc.integrated, reduction = "umap", group.by = "branch.anno")
options(repr.plot.width=16, repr.plot.height=12)
gl <- lapply(list(p1, p2, p3, p4), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}
sessionInfo()
R version 4.2.2 (2022-10-31) Platform: x86_64-conda-linux-gnu (64-bit) Running under: AlmaLinux 9.3 (Shamrock Pampas Cat) Matrix products: default BLAS/LAPACK: /hpc/group/pbenfeylab/ch416/miniconda3/envs/seu4/lib/libopenblasp-r0.3.21.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] grid stats4 stats graphics grDevices utils datasets [8] methods base other attached packages: [1] forcats_0.5.2 stringr_1.5.1 [3] dplyr_1.1.3 purrr_1.0.2 [5] readr_2.1.3 tidyr_1.3.0 [7] tibble_3.2.1 tidyverse_1.3.2 [9] gplots_3.1.3 RColorBrewer_1.1-3 [11] scater_1.26.1 scuttle_1.8.0 [13] scattermore_1.2 ggplot2_3.4.4 [15] cowplot_1.1.1 tricycle_1.6.0 [17] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0 [19] Biobase_2.58.0 GenomicRanges_1.50.0 [21] GenomeInfoDb_1.34.8 IRanges_2.32.0 [23] S4Vectors_0.36.0 BiocGenerics_0.44.0 [25] MatrixGenerics_1.10.0 matrixStats_1.1.0 [27] SeuratObject_4.1.3 Seurat_4.1.1.9001 loaded via a namespace (and not attached): [1] utf8_1.2.4 spatstat.explore_3.2-5 [3] reticulate_1.34.0 tidyselect_1.2.0 [5] RSQLite_2.3.1 AnnotationDbi_1.60.0 [7] htmlwidgets_1.6.2 BiocParallel_1.32.5 [9] Rtsne_0.16 munsell_0.5.0 [11] ScaledMatrix_1.6.0 codetools_0.2-19 [13] ica_1.0-3 pbdZMQ_0.3-8 [15] future_1.33.0 miniUI_0.1.1.1 [17] withr_2.5.2 spatstat.random_3.2-1 [19] colorspace_2.1-0 progressr_0.14.0 [21] uuid_1.1-0 ROCR_1.0-11 [23] tensor_1.5 listenv_0.9.0 [25] repr_1.1.4 GenomeInfoDbData_1.2.9 [27] polyclip_1.10-6 bit64_4.0.5 [29] parallelly_1.36.0 vctrs_0.6.4 [31] generics_0.1.3 circular_0.4-95 [33] timechange_0.1.1 R6_2.5.1 [35] ggbeeswarm_0.7.1 rsvd_1.0.5 [37] bitops_1.0-7 spatstat.utils_3.0-4 [39] cachem_1.0.8 DelayedArray_0.24.0 [41] assertthat_0.2.1 promises_1.2.1 [43] scales_1.2.1 googlesheets4_1.0.1 [45] beeswarm_0.4.0 gtable_0.3.4 [47] beachmat_2.14.0 globals_0.16.2 [49] goftest_1.2-3 rlang_1.1.2 [51] splines_4.2.2 lazyeval_0.2.2 [53] gargle_1.2.1 broom_1.0.2 [55] spatstat.geom_3.2-7 modelr_0.1.10 [57] reshape2_1.4.4 abind_1.4-5 [59] backports_1.4.1 httpuv_1.6.12 [61] tools_4.2.2 ellipsis_0.3.2 [63] ggridges_0.5.4 Rcpp_1.0.11 [65] plyr_1.8.9 base64enc_0.1-3 [67] sparseMatrixStats_1.10.0 zlibbioc_1.44.0 [69] RCurl_1.98-1.6 deldir_1.0-9 [71] pbapply_1.7-2 viridis_0.6.4 [73] zoo_1.8-12 haven_2.5.1 [75] ggrepel_0.9.4 cluster_2.1.4 [77] fs_1.6.3 magrittr_2.0.3 [79] data.table_1.14.8 RSpectra_0.16-1 [81] reprex_2.0.2 lmtest_0.9-40 [83] RANN_2.6.1 googledrive_2.0.0 [85] mvtnorm_1.1-3 ggnewscale_0.4.8 [87] fitdistrplus_1.1-11 hms_1.1.2 [89] patchwork_1.1.3 mime_0.12 [91] evaluate_0.23 xtable_1.8-4 [93] readxl_1.4.1 fastDummies_1.7.3 [95] gridExtra_2.3 compiler_4.2.2 [97] KernSmooth_2.23-20 crayon_1.5.2 [99] htmltools_0.5.7 tzdb_0.3.0 [101] later_1.3.1 lubridate_1.9.0 [103] DBI_1.1.3 dbplyr_2.2.1 [105] MASS_7.3-58.3 boot_1.3-28.1 [107] Matrix_1.5-4 cli_3.6.1 [109] parallel_4.2.2 igraph_1.5.1 [111] pkgconfig_2.0.3 sp_2.1-1 [113] IRdisplay_1.1 plotly_4.10.3 [115] spatstat.sparse_3.0-3 xml2_1.3.3 [117] vipor_0.4.5 XVector_0.38.0 [119] rvest_1.0.3 digest_0.6.33 [121] sctransform_0.4.1 RcppAnnoy_0.0.21 [123] spatstat.data_3.0-3 Biostrings_2.66.0 [125] cellranger_1.1.0 leiden_0.4.3 [127] uwot_0.1.16 DelayedMatrixStats_1.20.0 [129] shiny_1.7.5.1 gtools_3.9.4 [131] lifecycle_1.0.4 nlme_3.1-162 [133] jsonlite_1.8.7 BiocNeighbors_1.16.0 [135] viridisLite_0.4.2 fansi_1.0.5 [137] pillar_1.9.0 lattice_0.21-8 [139] KEGGREST_1.38.0 fastmap_1.1.1 [141] httr_1.4.7 survival_3.4-0 [143] glue_1.6.2 png_0.1-8 [145] bit_4.0.5 stringi_1.8.1 [147] blob_1.2.4 RcppHNSW_0.5.0 [149] BiocSingular_1.14.0 caTools_1.18.2 [151] memoise_2.0.1 IRkernel_1.3.1.9000 [153] irlba_2.3.5.1 future.apply_1.11.0
cc.genes.goldy <- read.csv("./tricycle/Goldy_2021_CycB1-1-sorted_Enriched_Depleted_G2M_rm-proto-genes.csv")
cc.genes.zhang <- read.csv("./tricycle/Zhang_2021_Dev_cell_Core_cell_cycle_genes.csv")
cc.genes.amigo.go <- as.character(read.table("./tricycle/GO_cell_cycle_annotation_arabidopsis.txt", header=F)$V1)
cc.genes.f <- unique(c(cc.genes.zhang$Gene, cc.genes.goldy$gene, cc.genes.amigo.go))
length(cc.genes.f)
pp.genes <- as.character(read.table("./CW_data/escoring/Root_sc/Protoplasting_DEgene_FC2_list.txt", header=F)$V1)
table(cc.genes.zhang$Phase)
G1/G0 G2 M S 18 95 77 53
head(cc.genes.zhang)
| Gene | Phase | Symbol | Araport11 | Function | |
|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | |
| 1 | AT1G12845 | G1/G0 | NA | transmembrane protein;(source:Araport11) | NA |
| 2 | AT3G26960 | G1/G0 | NA | Pollen Ole e 1 allergen and extensin family protein;(source:Araport11) | NA |
| 3 | AT5G02260 | G1/G0 | ATEXP9 | expansin A9;(source:Araport11) | member of Alpha-Expansin Gene Family. Naming convention from the Expansin Working Group (Kende et al, 2004. Plant Mol Bio) |
| 4 | AT1G33930 | G1/G0 | NA | P-loop containing nucleoside triphosphate hydrolases superfamily protein;(source:Araport11) | NA |
| 5 | AT5G01445 | G1/G0 | NA | hypothetical protein;(source:Araport11) | NA |
| 6 | AT5G67260 | G1/G0 | CYCD3;2 | CYCLIN D3;(source:Araport11) | Encode CYCD3;2, a CYCD3 D-type cyclin. Important for determining cell number in developing lateral organs. Mediating cytokinin effects in apical growth and development. |
## Cell Cyle genes used
cc.gene.list <- as.list(cc.genes.f)
rc.integrated <- readRDS("./scRNA-seq/Integrated_Objects/rc.integrated_18S_WT_cell_cycle_seu4_245genes_Goldy_Zhang_AMIGO_20230825.rds")
rownames(rc.integrated)
length(rownames(rc.integrated))
write.csv(data.frame("GeneID"=rownames(rc.integrated)),"./tradeseq/245_cell_cycle_related_genes_for_reference.csv")
## Prepare New Ref (All cells)
corrected.m <- rc.integrated@assays$integrated@data
# run PCA on the batch effects corrected matrix and get the rotaions scores for the top 2 PCs
pca.m <- scater::calculatePCA(corrected.m, ntop = 500)
new.ref <- attr(pca.m, 'rotation')[, seq_len(2)]
saveRDS(new.ref, "./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")
unique(rc.integrated$orig.ident)
#suppressMessages(rc.integrated <- RunUMAP(rc.integrated, reduction = "pca", dims = 1:8))
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "seurat_clusters", label = TRUE)+theme(legend.position="none")
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "seurat_clusters", label = TRUE)+theme(legend.position="none")
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "orig.ident", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "orig.ident", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Xylem Pole Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li <- factor(rc.integrated$celltype.anno.Li, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li),order))]
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li", label = FALSE, cols=color, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "celltype.anno.Li", label = FALSE, cols=color, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "time.anno.Li", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "consensus.time.group", label = FALSE, pt.size=1, order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
## S-phase cluster
sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==0)]$celltype.anno.Li))
Quiescent Center Phloem Phloem Pole Pericycle
38 275 323
Columella Xylem Cortex
331 795 1106
Procambium Xylem Pole Pericycle Endodermis
1426 1439 1592
Atrichoblast Trichoblast Lateral Root Cap
1646 2546 2879
## S-phase cluster
sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==0)]$time.anno.Li))
Distal Lateral Root Cap Distal Columella Proximal Columella
136 161 170
Maturation Transition Domain Proximal Lateral Root Cap
763 1900 2743
Proliferation Domain Elongation
3610 4913
## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$celltype.anno.Li)))
Lateral Root Cap Atrichoblast Endodermis
2981 1013 790
Trichoblast Cortex Xylem
550 429 284
Columella Xylem Pole Pericycle Procambium
275 252 249
Phloem Phloem Pole Pericycle Quiescent Center
233 100 8
## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$time.anno.Li)))
Proximal Lateral Root Cap Proliferation Domain Elongation
2863 1840 1324
Maturation Transition Domain Distal Columella
533 211 185
Distal Lateral Root Cap Proximal Columella
118 90
## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$time.celltype.anno.Li)))
Proximal Lateral Root Cap
2863
Proliferation Domain_Atrichoblast
828
Proliferation Domain_Endodermis
509
Proliferation Domain_Cortex
264
Proliferation Domain_Trichoblast
210
Elongation_Xylem Pole Pericycle
196
Elongation_Endodermis
195
Distal Columella
185
Elongation_Procambium
176
Elongation_Metaphloem & Companion Cell
153
Transition Domain_Trichoblast
144
Elongation_Protoxylem
128
Maturation_Protoxylem
123
Elongation_Cortex
118
Distal Lateral Root Cap
118
Elongation_Trichoblast
114
Elongation_Atrichoblast
99
Proximal Columella
90
Elongation_Phloem Pole Pericycle
90
Maturation_Endodermis
85
Maturation_Trichoblast
82
Maturation_Procambium
71
Maturation_Atrichoblast
57
Maturation_Xylem Pole Pericycle
41
Maturation_Cortex
33
Transition Domain_Atrichoblast
29
Elongation_Protophloem
29
Maturation_Metaphloem & Companion Cell
28
Elongation_Metaxylem
26
Transition Domain_Metaphloem & Companion Cell
14
Transition Domain_Cortex
14
Proliferation Domain_Xylem Pole Pericycle
12
Proliferation Domain_Quiescent Center
8
Maturation_Phloem Pole Pericycle
8
Proliferation Domain_Metaxylem
5
Maturation_Protophloem
5
Transition Domain_Xylem Pole Pericycle
3
Transition Domain_Protophloem
2
Transition Domain_Phloem Pole Pericycle
2
Proliferation Domain_Metaphloem & Companion Cell
2
Transition Domain_Procambium
1
Transition Domain_Metaxylem
1
Transition Domain_Endodermis
1
Proliferation Domain_Protoxylem
1
Proliferation Domain_Procambium
1
pc1r <- as.data.frame(rc.integrated[['pca']][,1]) %>% arrange(desc(PC_1)) %>% rownames(.)
pc1r[1:9]
## PC1 positive
FeaturePlot(rc.integrated, reduction="pca", features = c(pc1r[1:9]), cols = c("#EEEEEE","#FF4040"), order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
## PC1 positive
FeaturePlot(rc.integrated, reduction="umap", features = c(pc1r[1:9]), cols = c("#EEEEEE","#FF4040"), order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
pc1r[1:20]
## PC1 negative
FeaturePlot(rc.integrated, reduction="pca", features = c(rev(pc1r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
pc2r <- as.data.frame(rc.integrated[['pca']][,2]) %>% arrange(desc(PC_2)) %>% rownames(.)
pc2r[1:20]
rev(pc2r)[1:20]
## PC2 positive
FeaturePlot(rc.integrated, reduction="pca", features = c(pc2r[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
## PC2 negative
FeaturePlot(rc.integrated, reduction="pca", features = c(rev(pc2r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
## PC2 negative
FeaturePlot(rc.integrated, reduction="umap", features = c(rev(pc2r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
pc3r <- as.data.frame(rc.integrated[['pca']][,3]) %>% arrange(desc(PC_3)) %>% rownames(.)
pc3r[1:20]
rev(pc3r)[1:20]
zscore <- function(x){(x-mean(x))/sd(x)}
cc.genes.zhang <- cc.genes.zhang[match(intersect(rownames(rc.integrated),cc.genes.zhang$Gene),cc.genes.zhang$Gene),] %>% arrange(Phase)
head(cc.genes.zhang)
| Gene | Phase | Symbol | Araport11 | Function | |
|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | |
| 1 | AT1G16630 | G2 | NA | transmembrane protein;(source:Araport11) | NA |
| 2 | AT4G05190 | G2 | ATK5 | kinesin 5;(source:Araport11) | ATK5 encodes a kinesin protein involved in microtubule spindle morphogenesis. It acts as a minus-end directed motor as well as a plus-end tracking protein (+TIP). Localizes to mitotic spindle midzones and regions rich in growing plus-ends within phragmoplasts. |
| 3 | AT1G20610 | G2 | CYCB2;3 | Cyclin B2;(source:Araport11) | NA |
| 4 | AT5G55520 | G2 | NA | kinesin-like protein;(source:Araport11) | NA |
| 5 | AT4G01730 | G2 | NA | DHHC-type zinc finger family protein;(source:Araport11) | NA |
| 6 | AT1G44110 | G2 | CYCA1;1 | Cyclin A1;(source:Araport11) | NA |
head(cc.genes.zhang %>% filter(Phase == "S"))
| Gene | Phase | Symbol | Araport11 | Function | |
|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | |
| 1 | AT5G10400 | S | H3.1 | Histone superfamily protein;(source:Araport11) | NA |
| 2 | AT5G22880 | S | H2B | histone B2;(source:Araport11) | Encodes a histone 2B (H2B) protein. This protein can be ubiquitinated in planta, and this modification depends on the HUB1 and HUB2 E3 ubiquitin ligases. |
| 3 | AT5G59870 | S | h2a.w.6 | histone H2A 6;(source:Araport11) | Encodes HTA6, a histone H2A protein. |
| 4 | AT5G59690 | S | NA | Histone superfamily protein;(source:Araport11) | NA |
| 5 | AT5G65360 | S | H3.1 | Histone superfamily protein;(source:Araport11) | NA |
| 6 | AT1G09200 | S | H3.1 | Histone superfamily protein;(source:Araport11) | NA |
msc <- c()
suppressWarnings(
for (i in as.character(unique(cc.genes.zhang$Phase))){
if (length(cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene)>1){
msc <- cbind(msc, as.numeric(apply(apply(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene,], 1, zscore), 1, mean)))
} else {
msc <- cbind(msc, as.numeric(zscore(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene,])))
}
}
)
colnames(msc) <- as.character(unique(cc.genes.zhang$Phase))
rownames(msc) <- colnames(rc.integrated)
head(msc)
| G2 | M | S | |
|---|---|---|---|
| AAACCTGAGCGCCTCA_1 | -0.028233902 | 0.02767077 | -0.2789407 |
| AAACCTGCACTTAACG_1 | -0.085037611 | -0.09186232 | -0.3035347 |
| AAACCTGCAGGACGTA_1 | -0.128067505 | -0.03061144 | -0.2891731 |
| AAACCTGGTATCACCA_1 | -0.005055825 | -0.07676504 | -0.2121654 |
| AAACCTGTCAGGTAAA_1 | -0.082317857 | -0.09259665 | -0.1922454 |
| AAACGGGCAATAGAGT_1 | -0.154851348 | -0.11182951 | -0.3160343 |
#rc.integrated$G1_sig <- msc[,1]
rc.integrated$G2_sig <- msc[,1]
rc.integrated$M_sig <- msc[,2]
rc.integrated$S_sig <- msc[,3]
#FeaturePlot(rc.integrated, reduction="umap", features = "G1_sig", pt.size=1,order=TRUE)
FeaturePlot(rc.integrated, reduction="umap", features = "S_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="umap", features = "G2_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="umap", features = "M_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
#FeaturePlot(rc.integrated, reduction="pca", features = "G1_sig", pt.size=1,order=TRUE)
FeaturePlot(rc.integrated, reduction="pca", features = "S_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="pca", features = "G2_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="pca", features = "M_sig", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
head(cc.genes.goldy)
| gene | phase | |
|---|---|---|
| <chr> | <chr> | |
| 1 | AT2G23320 | Enriched_G2M |
| 2 | AT4G17340 | Enriched_G2M |
| 3 | AT5G27420 | Enriched_G2M |
| 4 | AT5G54290 | Enriched_G2M |
| 5 | AT5G08240 | Enriched_G2M |
| 6 | AT5G06980 | Enriched_G2M |
cc.genes.goldy <- cc.genes.goldy[match(intersect(rownames(rc.integrated),cc.genes.goldy$gene),cc.genes.goldy$gene),] %>% arrange(phase)
msc <- c()
suppressWarnings(
for (i in as.character(unique(cc.genes.goldy$phase))){
if (length(cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene)>1){
msc <- cbind(msc, as.numeric(apply(apply(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene,], 1, zscore), 1, mean)))
} else {
msc <- cbind(msc, as.numeric(zscore(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene,])))
}
}
)
colnames(msc) <- as.character(unique(cc.genes.goldy$phase))
rownames(msc) <- colnames(rc.integrated)
head(msc)
| Depleted_G2M | Enriched_G2M | |
|---|---|---|
| AAACCTGAGCGCCTCA_1 | -0.055274931 | 0.05566898 |
| AAACCTGCACTTAACG_1 | 0.002564458 | -0.10651000 |
| AAACCTGCAGGACGTA_1 | -0.072215184 | -0.10086421 |
| AAACCTGGTATCACCA_1 | -0.161443430 | 0.01979224 |
| AAACCTGTCAGGTAAA_1 | 0.028097562 | -0.05966497 |
| AAACGGGCAATAGAGT_1 | 0.092315593 | -0.06260356 |
rc.integrated$Depleted_G2M <- msc[,1]
rc.integrated$Enriched_G2M <- msc[,2]
FeaturePlot(rc.integrated, reduction="umap", features = "Depleted_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="umap", features = "Enriched_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="pca", features = "Depleted_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
FeaturePlot(rc.integrated, reduction="pca", features = "Enriched_G2M", pt.size=1,order=TRUE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
rc.integrated$G2M_clust <- rep("others", ncol(rc.integrated))
rc.integrated$G2M_clust[which(rc.integrated$seurat_clusters==4)] = "G2M"
DimPlot(rc.integrated, reduction="umap", group = "G2M_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
DimPlot(rc.integrated, reduction="pca", group = "G2M_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
rc.integrated$S_clust <- rep("others", ncol(rc.integrated))
rc.integrated$S_clust[which(rc.integrated$seurat_clusters==0)] = "S"
DimPlot(rc.integrated, reduction="umap", group = "S_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
DimPlot(rc.integrated, reduction="pca", group = "S_clust", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
new.ref <- readRDS("./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")
sce <- as.SingleCellExperiment(rc.integrated)
sce$geno <- rep("WT",length(sce$orig.ident))
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
options(repr.plot.width=9, repr.plot.height=4)
print(plot_ccposition_den(sce$tricyclePosition,
sce$geno, 'genotype',
bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted",
color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted",
color = "red", size=1.5))
The number of projection genes found in the new data is 245.
Warning message:
“Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.”
grDevices::cairo_pdf("./pdfs/Density_Plot_CC.pdf",width=9, height=4)
print(plot_ccposition_den(sce$tricyclePosition,
sce$geno, 'genotype',
bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted",
color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted",
color = "red", size=1.5))
dev.off()
p <- plot_emb_circle_scale(sce, dimred = 1,
point.size = 5, point.alpha = 0.9) +
theme_bw(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
options(repr.plot.width=8, repr.plot.height=4)
plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.6))
grDevices::cairo_pdf("./pdfs/PCA_CC.pdf",width=8, height=4)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.6)))
dev.off()
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
options(repr.plot.width=9, repr.plot.height=9)
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted",
color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted",
color = "red", size=1.5))
The number of projection genes found in the new data is 245.
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
rc.integrated$tricyclePosition <- sce$tricyclePosition
rc.integrated$tricycleCCStage <- tricycleCCStage
saveRDS(rc.integrated,"./scRNA-seq/Integrated_Objects/rc.integrated_18S_WT_cell_cycle_seu4_245genes_Goldy_Zhang_AMIGO_20230825.rds")
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 1, point.alpha = 0.9) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
DimPlot(rc.integrated, reduction="umap", group = "tricycleCCStage", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
DimPlot(rc.integrated, reduction="pca", group = "tricycleCCStage", pt.size=1,order=FALSE)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
rc.integrated$CCclust <- rep("others", ncol(rc.integrated))
rc.integrated$CCclust[which(rc.integrated$seurat_clusters==0)] = "S"
rc.integrated$CCclust[which(rc.integrated$seurat_clusters==4)] = "G2M"
table(rc.integrated$tricycleCCStage, rc.integrated$CCclust)
G2M others S
G1/G0 920 154117 1510
G2M 6101 4126 117
S 143 20007 12769
scater::plotReducedDim(sce, dimred = "tricycleEmbedding") +
labs(x = "Projected PC1", y = "Projected PC2") +
ggtitle(sprintf("Projected cell cycle space (n=%d)",
ncol(sce))) +
theme_bw(base_size = 14)
p <- plot_emb_circle_scale(sce, dimred = "tricycleEmbedding",
point.size = 3.5, point.alpha = 0.9) +
theme_bw(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.4))
## Total RNA content
#fit.l <- fit_periodic_loess(sce$tricyclePosition,
# log2(sce$nCount_RNA),
# plot = TRUE, span=0.3,
# x_lab = "Cell cycle position \u03b8", y_lab = "log2(Total UMIs)",
# fig.title = paste0("Expression of key gene along \u03b8 (n=",
# ncol(sce), ")"))
#fit <- fit.l$pred.df
#options(repr.plot.width=8, repr.plot.height=6)
#ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "log2(total UMIs)", title = "") +
# scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0,
# pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0,
# 0.5, 1, 1.5, 2), "π"), name = "θ")+
# theme_bw(base_size = 14)+
# theme(legend.title=element_blank())
table(rc.integrated$time.anno.Li)
Distal Columella Distal Lateral Root Cap Elongation
6526 4904 83690
Maturation Proliferation Domain Proximal Columella
21054 21913 3063
Proximal Lateral Root Cap Transition Domain
27495 31165
sub <- subset(rc.integrated, cells=colnames(rc.integrated)[which(rc.integrated$time.anno.Li=="Proliferation Domain"|
rc.integrated$time.anno.Li=="Proximal Columella"|
rc.integrated$time.anno.Li=="Proximal Lateral Root Cap")])
sub
An object of class Seurat 57118 features across 52471 samples within 3 assays Active assay: integrated (245 features, 245 variable features) 2 layers present: data, scale.data 2 other assays present: RNA, SCT 3 dimensional reductions calculated: pca, umap, umap_2D
table(sub$time.anno.Li)
Proliferation Domain Proximal Columella Proximal Lateral Root Cap
21913 3063 27495
# Run PCA
sub <- RunPCA(sub, npcs = 10, verbose = FALSE, approx = TRUE)
# Find nearest neighbors
#rc.integrated <- FindNeighbors(rc.integrated, reduction = "pca",dims = 1:50)
Warning message: “Number of dimensions changing from 50 to 10” Warning message: “The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' This message will be shown once per session” 20:27:45 UMAP embedding parameters a = 0.9922 b = 1.112 Found more than one class "dist" in cache; using the first, from namespace 'spam' Also defined by ‘BiocGenerics’ 20:27:45 Read 52471 rows and found 5 numeric columns 20:27:45 Using Annoy for neighbor search, n_neighbors = 30 Found more than one class "dist" in cache; using the first, from namespace 'spam' Also defined by ‘BiocGenerics’ 20:27:45 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 20:27:57 Writing NN index file to temp file /tmp/Rtmpu0VjyG/file317db91a27d9c1 20:27:57 Searching Annoy index using 1 thread, search_k = 3000 20:28:40 Annoy recall = 100% 20:28:42 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 20:28:44 5 smooth knn distance failures 20:28:51 Initializing from normalized Laplacian + noise (using RSpectra) 20:28:59 Commencing optimization for 200 epochs, with 1889004 positive edges 20:29:42 Optimization finished Warning message: “Key ‘UMAP_’ taken, using ‘umap_’ instead”
# Run 2D UMAP
sub <- RunUMAP(sub, reduction = "pca", dims = 1:3, n.neighbors = 30, min.dist=0.3)
20:40:15 UMAP embedding parameters a = 0.9922 b = 1.112 Found more than one class "dist" in cache; using the first, from namespace 'spam' Also defined by ‘BiocGenerics’ 20:40:15 Read 52471 rows and found 3 numeric columns 20:40:15 Using Annoy for neighbor search, n_neighbors = 30 Found more than one class "dist" in cache; using the first, from namespace 'spam' Also defined by ‘BiocGenerics’ 20:40:15 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 20:40:28 Writing NN index file to temp file /tmp/Rtmpu0VjyG/file317db92623cbe 20:40:28 Searching Annoy index using 1 thread, search_k = 3000 20:40:57 Annoy recall = 100% 20:40:59 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 20:41:00 18 smooth knn distance failures 20:41:06 Initializing from normalized Laplacian + noise (using RSpectra) 20:41:40 Commencing optimization for 200 epochs, with 1402752 positive edges 20:42:13 Optimization finished Warning message: “Key ‘UMAP_’ taken, using ‘umap_’ instead”
my_cols = brewer.pal(8,"Set2")
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(sub, reduction = "umap", group.by = "time.anno.Li",cols=alpha(my_cols,0.8), label = FALSE, pt.size=1)
sce <- as.SingleCellExperiment(sub)
# dims = 3
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 2, point.alpha = 0.8) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
# dims = 4
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 2, point.alpha = 1) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
# dims = 5
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 2, point.alpha = 1) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
# dims = 5
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 2, point.alpha = 1) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
run_tricycle <- function(sample){
seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
sce <- as.SingleCellExperiment(seu)
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8") +
theme_bw(base_size = 14))
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
## Schwabe Stage
#sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
## Table
#seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
#print(table(seu$time.anno.Li, sce$tricycleCCStage))
## Prepare seu
#seu <- prep(seu)
#plot_anno(seu)
## Sce to seu
seu$tricyclePosition <- sce$tricyclePosition/pi
seu$tricycleCCStage <- sce$tricycleCCStage
#seu$SchwabeCCStage <- sce$CCStage
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 1, point.alpha = 0.9) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
## Write annotations as csv
write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
#saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}
sample.list <- c("sc_130","sc_131","sc_132","sc_133","sc_134","sc_135","sc_136","sc_137")
for (i in 1:length(sample.list)){
run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 230. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ” The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 237. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
sample.list <- c("sc_5")
for (i in 1:length(sample.list)){
run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
sample.list <- c("sc_1","sc_2","sc_111","sc_112","sc_113","sc_114")
for (i in 1:length(sample.list)){
run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ” The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 243. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
sample.list <- c("sc_151","sc_152","sc_153","sc_154")
for (i in 1:length(sample.list)){
run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ” The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
sample.list <- c("pp1","col0","tnw2","sc_1","sc_9_at","sc_30", "sc_31", "sc_37", "sc_40", "sc_51")
for (i in 1:length(sample.list)){
run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ” The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 243. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 244. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 241. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
run_tricycle <- function(sample){
seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
sce <- as.SingleCellExperiment(seu)
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.2
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- (sce$tricyclePosition-2)*-1
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8") +
theme_bw(base_size = 14))
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
## Schwabe Stage
#sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
## Table
#seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
#print(table(seu$time.anno.Li, sce$tricycleCCStage))
## Prepare seu
#seu <- prep(seu)
#plot_anno(seu)
## Sce to seu
seu$tricyclePosition <- sce$tricyclePosition/pi
seu$tricycleCCStage <- sce$tricycleCCStage
#seu$SchwabeCCStage <- sce$CCStage
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 1, point.alpha = 0.9) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
## Write annotations as csv
write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
#saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}
sample.list <- c("dc1","dc2","tnw1","briT","briTR","sc_10_at","sc_11","sc_12","sc_20","sc_43","sc_44","sc_45","sc_46","sc_47","sc_48","sc_49","sc_50","sc_71",
"sc_122","sc_123","sc_124","sc_125","sc_126","sc_127","sc_128","sc_129","sc_130","sc_131","sc_132","sc_133","sc_134","sc_135","sc_136","sc_137",
"sc_155","sc_156","sc_157","sc_158","sc_159","sc_160","sc_161","sc_162","sc_163","sc_164","sc_165","sc_166","sc_167","sc_168","sc_169",
"sc_170","sc_171","sc_172","sc_173","sc_175","sc_176","sc_177","sc_178","sc_179","sc_180","sc_186","sc_187","sc_189","sc_190","sc_191",
"sc_215","sc_216","sc_217","sc_218","sc_220","sc_221","sc_222","sc_223","sc_224","sc_225","sc_226","sc_227","sc_228","sc_229","sc_235")
for (i in 1:length(sample.list)){
run_tricycle(sample.list[i])
}
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ” The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 240. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 228. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 244. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
The number of projection genes found in the new data is 245. Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/","dc1","_COPILOT.rds"))
sce <- as.SingleCellExperiment(seu)
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
The number of projection genes found in the new data is 245.
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8") +
theme_bw(base_size = 14))
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
## Schwabe Stage
#sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
## Table
#seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
#print(table(seu$time.anno.Li, sce$tricycleCCStage))
## Prepare seu
#seu <- prep(seu)
#plot_anno(seu)
## Sce to seu
seu$tricyclePosition <- sce$tricyclePosition/pi
seu$tricycleCCStage <- sce$tricycleCCStage
#seu$SchwabeCCStage <- sce$CCStage
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 1, point.alpha = 0.9) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
run_tricycle <- function(sample){
seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
sce <- as.SingleCellExperiment(seu)
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8") +
theme_bw(base_size = 14))
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
## Schwabe Stage
#sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
## Table
#seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
#print(table(seu$time.anno.Li, sce$tricycleCCStage))
## Prepare seu
#seu <- prep(seu)
#plot_anno(seu)
## Sce to seu
seu$tricyclePosition <- sce$tricyclePosition/pi
seu$tricycleCCStage <- sce$tricycleCCStage
#seu$SchwabeCCStage <- sce$CCStage
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 1, point.alpha = 0.9) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
## Write annotations as csv
write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
#saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}
mtx <- read.table('./scRNA-seq/Published_Data/GSE140778_20190115_data.txt')
meta <- read.csv('./scRNA-seq/Published_Data/GSE140778_20190115_meta_data.csv')
sum(!is.na(meta$Cluster_PSE))
head(mtx)
| SLX_12100_i701_i502 | SLX_12100_i701_i503 | SLX_12100_i701_i504 | SLX_12100_i701_i505 | SLX_12100_i701_i506 | SLX_12100_i701_i507 | SLX_12100_i701_i508 | SLX_12100_i701_i517 | SLX_12100_i702_i502 | SLX_12100_i702_i503 | ... | SLX_17313_i728_i521 | SLX_17313_i728_i522 | SLX_17313_i729_i513 | SLX_17313_i729_i515 | SLX_17313_i729_i516 | SLX_17313_i729_i517 | SLX_17313_i729_i518 | SLX_17313_i729_i520 | SLX_17313_i729_i521 | SLX_17313_i729_i522 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | ... | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| AT1G01010 | 138.832 | 0.000 | 0 | 1.37359 | 0.665616 | 23.122100 | 0.0000 | 0 | 0.00000 | 153.1940 | ... | 45.11770 | 0.00000 | 0 | 13.7543 | 58.2903 | 0.00 | 0.00000 | 0 | 0.0000 | 0 |
| AT1G01020 | 277.664 | 146.593 | 0 | 33.42390 | 0.000000 | 0.000000 | 39.6054 | 0 | 0.00000 | 105.3210 | ... | 8.01155 | 0.00000 | 0 | 0.0000 | 0.0000 | 0.00 | 0.00000 | 0 | 0.0000 | 0 |
| AT1G01030 | 0.000 | 0.000 | 0 | 0.00000 | 7.321770 | 0.000000 | 0.0000 | 0 | 0.00000 | 0.0000 | ... | 0.00000 | 0.00000 | 0 | 0.0000 | 0.0000 | 0.00 | 0.00000 | 0 | 0.0000 | 0 |
| AT1G01040 | 0.000 | 0.000 | 0 | 0.00000 | 0.332834 | 0.379051 | 22.9967 | 0 | 0.00000 | 4.7873 | ... | 0.00000 | 2.43814 | 0 | 0.0000 | 26.1001 | 0.00 | 0.00000 | 0 | 0.0000 | 0 |
| AT1G01046 | 0.000 | 0.000 | 0 | 0.00000 | 0.000000 | 0.000000 | 0.0000 | 0 | 0.00000 | 0.0000 | ... | 0.00000 | 0.00000 | 0 | 0.0000 | 0.0000 | 0.00 | 0.00000 | 0 | 0.0000 | 0 |
| AT1G01050 | 138.832 | 146.593 | 0 | 72.80000 | 105.167000 | 137.974000 | 136.7020 | 0 | 1.94657 | 4.7873 | ... | 6.32491 | 104.80900 | 0 | 103.1570 | 29.5801 | 329.67 | 6.51429 | 0 | 10.5481 | 0 |
rownames(meta) <- meta$ID
seu <- CreateSeuratObject(mtx, meta.data=meta)
seu
An object of class Seurat 32824 features across 1375 samples within 1 assay Active assay: RNA (32824 features, 0 variable features)
head(seu@meta.data)
| orig.ident | nCount_RNA | nFeature_RNA | ID | sample | reads_total | reads_unmapped | reads_unique | reads_multi | percent_unique | percent_multi | percent_mapped | lib_size_avr | Cluster_all_cells | Cluster_PSE | Pseudotime_PSE | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <fct> | <dbl> | <int> | <chr> | <chr> | <int> | <int> | <int> | <int> | <dbl> | <dbl> | <dbl> | <int> | <int> | <int> | <dbl> | |
| SLX_12100_i701_i502 | SLX | 734281.9 | 2860 | SLX_12100_i701_i502 | P1 | 847809 | 803895 | 38147 | 5767 | 4.499480 | 0.6802240 | 5.179704 | 678 | NA | NA | NA |
| SLX_12100_i701_i503 | SLX | 913422.6 | 3518 | SLX_12100_i701_i503 | P1 | 1431659 | 1399319 | 27698 | 4642 | 1.934679 | 0.3242392 | 2.258918 | 678 | NA | NA | NA |
| SLX_12100_i701_i504 | SLX | 827476.9 | 2904 | SLX_12100_i701_i504 | P1 | 313281 | 300049 | 9440 | 3792 | 3.013269 | 1.2104149 | 4.223684 | 678 | NA | NA | NA |
| SLX_12100_i701_i505 | SLX | 600775.5 | 10725 | SLX_12100_i701_i505 | P1 | 3328241 | 992332 | 1515355 | 820554 | 45.530206 | 24.6542844 | 70.184491 | 678 | 10 | 4 | 1.3333992 |
| SLX_12100_i701_i506 | SLX | 922948.5 | 11363 | SLX_12100_i701_i506 | P1 | 7051643 | 1901022 | 3201903 | 1948718 | 45.406482 | 27.6349498 | 73.041432 | 678 | 14 | 4 | 0.8876946 |
| SLX_12100_i701_i507 | SLX | 894200.0 | 10634 | SLX_12100_i701_i507 | P1 | 6906310 | 1881499 | 2873579 | 2151232 | 41.608022 | 31.1487900 | 72.756812 | 678 | 7 | NA | NA |
seu <- NormalizeData(seu)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = nrow(seu))
all.genes <- rownames(seu)
seu <- ScaleData(seu, features = all.genes)
seu <- RunPCA(seu,features = VariableFeatures(object = seu))
Centering and scaling data matrix PC_ 1 Positive: AT5G03240, AT3G53990, AT4G05320, AT1G70670, AT4G05050, AT1G75280, AT4G39090, AT3G22400, AT4G21980, AT1G29520 AT3G03200, AT1G77700, AT5G48060, AT1G54400, AT3G06420, AT3G58720, AT4G33490, AT2G46630, AT1G29380, AT1G08160 AT1G79480, AT5G07220, AT5G49320, AT5G18970, AT5G03510, AT2G41200, AT2G42820, AT1G76970, AT1G63310, AT1G61760 Negative: AT3G04400, AT3G24830, AT3G25520, AT5G20290, AT5G02960, AT1G08360, AT5G16130, AT1G18540, AT3G02080, AT3G49910 AT4G27090, AT2G34480, AT5G15200, AT2G47610, AT5G10360, AT4G00100, AT4G34670, AT5G23740, AT2G39460, AT3G04920 AT3G53020, AT1G69620, AT2G18020, AT2G27530, AT1G67430, AT3G14600, AT3G11940, AT3G09200, AT2G19730, AT4G18100 PC_ 2 Positive: AT5G42050, AT1G20450, AT1G19180, AT1G20440, AT2G47730, AT1G78380, AT4G30960, AT3G10985, AT5G10830, AT5G11670 AT5G45110, AT1G78660, AT5G63790, AT1G04250, AT2G29440, AT5G08790, AT4G16380, AT4G34150, AT5G64260, AT5G19110 AT2G24850, AT4G35090, AT4G19200, AT3G52470, AT5G06320, AT4G21990, AT1G55920, AT1G78080, AT1G76680, AT4G36040 Negative: AT5G62960, AT1G80940, AT2G38550, AT1G05610, AT2G22730, AT1G31880, AT5G03530, AT1G55760, AT2G18380, AT1G62045 AT1G79430, AT5G40580, AT1G63310, AT1G68000, AT1G61760, AT3G49360, AT4G23900, AT1G69970, AT2G20515, AT2G36400 AT4G17510, AT5G26940, AT1G14900, AT5G04060, AT1G06490, AT2G32280, AT5G17260, AT5G48060, AT4G24250, AT1G79480 PC_ 3 Positive: AT1G03220, AT1G48920, AT1G66580, AT5G52470, AT3G15000, AT5G20160, AT3G10690, AT1G74560, AT4G12600, AT1G56110 AT5G27120, AT4G26110, AT2G38810, AT5G25370, AT4G25630, AT1G68230, AT3G05060, AT3G07770, AT3G44750, AT5G22650 AT3G45970, AT1G52930, AT3G23990, AT3G58350, AT3G56070, AT2G19540, AT3G57150, AT4G15910, AT3G03950, AT1G26680 Negative: AT1G78040, AT5G60920, AT5G59880, AT2G45470, AT3G12610, AT5G39320, AT5G05170, AT5G49720, AT1G53500, AT1G75680 AT1G08200, AT4G12730, AT5G12250, AT2G20750, AT1G26850, AT1G70490, AT1G19835, AT4G31590, AT3G12710, AT1G10630 AT1G76160, AT4G34490, AT1G45688, AT4G32410, AT1G76670, AT5G59290, AT5G64740, AT5G47770, AT5G25100, AT3G07010 PC_ 4 Positive: AT1G10760, AT4G01610, AT1G54030, AT1G32860, AT2G34020, AT1G18210, AT2G29980, AT3G55760, AT5G48300, AT3G12490 AT4G09020, AT4G31290, AT1G21340, AT3G08690, AT1G56680, AT2G15490, AT2G32430, AT2G04025, AT4G00490, AT5G08590 AT2G28900, AT3G46440, AT5G18470, AT2G41800, AT4G34970, AT1G72490, AT1G13620, AT4G23060, AT3G18230, AT1G19850 Negative: AT4G14130, AT1G62480, AT4G21960, AT3G53420, AT2G02130, AT5G59090, AT1G52760, AT3G21770, AT2G43910, AT3G04730 AT4G29100, AT3G19450, AT3G23820, AT3G51950, AT5G56540, AT3G14310, AT4G11210, AT3G56240, AT2G30520, AT3G16330 AT5G24800, AT4G19410, AT2G31083, AT5G51780, AT2G44790, AT1G10682, AT3G16340, AT3G11930, AT5G25890, AT4G36710 PC_ 5 Positive: AT4G15960, AT1G70670, AT3G02470, AT4G18710, AT2G32150, AT2G39420, AT2G36650, AT3G14780, AT4G37730, AT5G16650 AT3G47160, AT4G35450, AT4G35790, AT2G38800, AT5G13740, AT3G11280, AT4G39090, AT2G19110, AT5G64572, AT1G75380 AT5G49320, AT1G80040, AT1G44170, AT5G18980, AT3G60680, AT3G63000, AT3G14350, AT4G21110, AT2G39080, AT2G21270 Negative: ATMG01390, ATMG00020, AT5G18610, AT2G07698, AT1G77690, AT5G23860, ATMG01190, AT1G04680, AT5G26270, AT3G02230 AT2G04780, ATMG01360, AT1G64640, AT3G54810, AT1G69690, AT2G14890, ATCG00500, AT2G20750, ATMG00660, AT2G47650 ATCG01180, ATCG00470, AT3G48185, AT2G21160, AT3G58120, ATCG00950, AT5G15650, ATMG00090, AT2G34300, AT1G21750
seu <- FindNeighbors(seu, dims = 1:10)
seu <- FindClusters(seu, resolution = 0.5)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 1375 Number of edges: 33789 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9051 Number of communities: 12 Elapsed time: 0 seconds
seu <- RunUMAP(seu, dims = 1:10)
Warning message: “The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' This message will be shown once per session” 10:49:41 UMAP embedding parameters a = 0.9922 b = 1.112 10:49:41 Read 1375 rows and found 10 numeric columns 10:49:41 Using Annoy for neighbor search, n_neighbors = 30 10:49:41 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 10:49:41 Writing NN index file to temp file /tmp/RtmpSzDebc/file16ca8b2394d636 10:49:41 Searching Annoy index using 1 thread, search_k = 3000 10:49:42 Annoy recall = 100% 10:49:43 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 10:49:46 Initializing from normalized Laplacian + noise (using RSpectra) 10:49:46 Commencing optimization for 500 epochs, with 49356 positive edges 10:49:49 Optimization finished
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(seu, reduction = "umap")
head(seu@meta.data)
| orig.ident | nCount_RNA | nFeature_RNA | ID | sample | reads_total | reads_unmapped | reads_unique | reads_multi | percent_unique | percent_multi | percent_mapped | lib_size_avr | Cluster_all_cells | Cluster_PSE | Pseudotime_PSE | RNA_snn_res.0.5 | seurat_clusters | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <fct> | <dbl> | <int> | <chr> | <chr> | <int> | <int> | <int> | <int> | <dbl> | <dbl> | <dbl> | <int> | <int> | <int> | <dbl> | <fct> | <fct> | |
| SLX_12100_i701_i502 | SLX | 734281.9 | 2860 | SLX_12100_i701_i502 | P1 | 847809 | 803895 | 38147 | 5767 | 4.499480 | 0.6802240 | 5.179704 | 678 | NA | NA | NA | 10 | 10 |
| SLX_12100_i701_i503 | SLX | 913422.6 | 3518 | SLX_12100_i701_i503 | P1 | 1431659 | 1399319 | 27698 | 4642 | 1.934679 | 0.3242392 | 2.258918 | 678 | NA | NA | NA | 4 | 4 |
| SLX_12100_i701_i504 | SLX | 827476.9 | 2904 | SLX_12100_i701_i504 | P1 | 313281 | 300049 | 9440 | 3792 | 3.013269 | 1.2104149 | 4.223684 | 678 | NA | NA | NA | 0 | 0 |
| SLX_12100_i701_i505 | SLX | 600775.5 | 10725 | SLX_12100_i701_i505 | P1 | 3328241 | 992332 | 1515355 | 820554 | 45.530206 | 24.6542844 | 70.184491 | 678 | 10 | 4 | 1.3333992 | 0 | 0 |
| SLX_12100_i701_i506 | SLX | 922948.5 | 11363 | SLX_12100_i701_i506 | P1 | 7051643 | 1901022 | 3201903 | 1948718 | 45.406482 | 27.6349498 | 73.041432 | 678 | 14 | 4 | 0.8876946 | 8 | 8 |
| SLX_12100_i701_i507 | SLX | 894200.0 | 10634 | SLX_12100_i701_i507 | P1 | 6906310 | 1881499 | 2873579 | 2151232 | 41.608022 | 31.1487900 | 72.756812 | 678 | 7 | NA | NA | 7 | 7 |
seu
An object of class Seurat 32824 features across 1375 samples within 1 assay Active assay: RNA (32824 features, 27218 variable features) 2 dimensional reductions calculated: pca, umap
#seu[["umap"]]@cell.embeddings[,1] <- seu[["umap"]]@cell.embeddings[,1]*-1
seu[["umap"]]@cell.embeddings[,2] <- seu[["umap"]]@cell.embeddings[,2]*-1
#u2 <- rc.integrated@reductions$umap@cell.embeddings[,1]
#u1 <- rc.integrated@reductions$umap@cell.embeddings[,2]
#rc.integrated@reductions$umap@cell.embeddings[,1] <- u1
#rc.integrated@reductions$umap@cell.embeddings[,2] <- u2
DimPlot(seu, reduction = "umap", group.by = "Cluster_PSE")
FeaturePlot(seu, reduction = "umap", features = "Pseudotime_PSE")
sce <- as.SingleCellExperiment(seu)
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
The number of projection genes found in the new data is 245.
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8") +
theme_bw(base_size = 14))
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
## Schwabe Stage
#sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
## Table
#seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
#print(table(seu$time.anno.Li, sce$tricycleCCStage))
## Prepare seu
#seu <- prep(seu)
#plot_anno(seu)
## Sce to seu
seu$tricyclePosition <- sce$tricyclePosition/pi
seu$tricycleCCStage <- sce$tricycleCCStage
#seu$SchwabeCCStage <- sce$CCStage
Warning message in brewer.pal(nlevels(color_var.v), "Set1"): “minimal value for n is 3, returning requested palette with 3 different levels ”
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 4, point.alpha = 1) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
DimPlot(seu, reduction = "umap", group.by = "tricycleCCStage")
seu$Pseudotime_PSE
seu$Pseudotime_PSE_rank <- rank(seu$Pseudotime_PSE)
FeaturePlot(seu, reduction = "umap", features = "Pseudotime_PSE_rank")
saveRDS(seu,"./tricycle/seu_Phloem_1375_cells.rds")
seu_PSE <- subset(seu, cells=colnames(seu)[!is.na(seu$Cluster_PSE)])
seu_PSE
An object of class Seurat 32824 features across 758 samples within 1 assay Active assay: RNA (32824 features, 27218 variable features) 2 dimensional reductions calculated: pca, umap
DimPlot(seu_PSE, reduction = "umap", group.by = "tricycleCCStage")
FeaturePlot(seu_PSE, reduction = "umap", features = "Pseudotime_PSE")
FeaturePlot(seu_PSE, reduction = "umap", features = "Pseudotime_PSE_rank")
saveRDS(seu_PSE,"./tricycle/seu_PSE_758_cells.rds")
## All cells (PSE, MSE, Procambium)
pltsg <- function(gene){
keygene.idx <- which(rownames(seu@assays$RNA@data) == gene)
fit.l <- fit_periodic_loess(seu$tricyclePosition*pi,
seu@assays$RNA@data[keygene.idx,],
plot = FALSE, span=0.3,
x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
fit <- fit.l$pred.df
options(repr.plot.width=8, repr.plot.height=6)
print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) +
scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0,
pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0,
0.5, 1, 1.5, 2), "π"), name = "θ")+
# scale_x_continuous(limits = c(0, 200), breaks = c(0,
# 50, 100, 150, 200), labels = paste0(c(0,
# 0.5, 1, 1.5, 2), "π"), name = "θ")+
theme_bw(base_size = 14)+
theme(legend.title=element_blank()))
}
## DWF4
pltsg('AT3G50660')
## OPS
pltsg('AT3G09070')
## CPD
pltsg('AT5G05690')
## Only PSE
pltsg <- function(gene){
keygene.idx <- which(rownames(seu_PSE@assays$RNA@data) == gene)
fit.l <- fit_periodic_loess(seu_PSE$tricyclePosition*pi,
seu_PSE@assays$RNA@data[keygene.idx,],
plot = FALSE, span=0.3,
x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
fit <- fit.l$pred.df
options(repr.plot.width=8, repr.plot.height=6)
print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) +
scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0,
pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0,
0.5, 1, 1.5, 2), "π"), name = "θ")+
# scale_x_continuous(limits = c(0, 200), breaks = c(0,
# 50, 100, 150, 200), labels = paste0(c(0,
# 0.5, 1, 1.5, 2), "π"), name = "θ")+
theme_bw(base_size = 14)+
theme(legend.title=element_blank()))
}
## DWF4
pltsg('AT3G50660')
## OPS
pltsg('AT3G09070')
## CPD
pltsg('AT5G05690')
seu_mer_PSE <- subset(seu_PSE, cells=colnames(seu_PSE)[which(seu_PSE$Pseudotime_PSE_rank <= 410)])
seu_mer_PSE
An object of class Seurat 32824 features across 410 samples within 1 assay Active assay: RNA (32824 features, 27218 variable features) 2 dimensional reductions calculated: pca, umap
saveRDS(seu_mer_PSE,"./tricycle/seu_mer_PSE_410_cells.rds")
## Only meristem PSE
pltsg <- function(gene){
keygene.idx <- which(rownames(seu_mer_PSE@assays$RNA@data) == gene)
fit.l <- fit_periodic_loess(seu_mer_PSE$tricyclePosition*pi,
seu_mer_PSE@assays$RNA@data[keygene.idx,],
plot = FALSE, span=0.3,
x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
fit <- fit.l$pred.df
options(repr.plot.width=8, repr.plot.height=6)
print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) +
scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0,
pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0,
0.5, 1, 1.5, 2), "π"), name = "θ")+
# scale_x_continuous(limits = c(0, 200), breaks = c(0,
# 50, 100, 150, 200), labels = paste0(c(0,
# 0.5, 1, 1.5, 2), "π"), name = "θ")+
theme_bw(base_size = 14)+
theme(legend.title=element_blank()))
}
## DWF4
pltsg('AT3G50660')
## OPS
pltsg('AT3G09070')
## CPD
pltsg('AT5G05690')
sce <- readRDS("./scRNA-seq/Published_Data/phloem_atlas_sce.rds")
meta <- read.csv('./scRNA-seq/Published_Data/phloem_atlas_cell_metadata.csv')
head(meta)
| cell_id | sample | barcode | total_counts | detected_genes | cluster | annotation | UMAP1 | UMAP2 | |
|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <int> | <int> | <int> | <chr> | <dbl> | <dbl> | |
| 1 | APL_AAACCCAAGTAACGTA-1 | APL | AAACCCAAGTAACGTA-1 | 35600 | 5210 | 9 | early/undetermined | 2.2210797 | -3.4197492 |
| 2 | APL_AAACCCACAACCAACT-1 | APL | AAACCCACAACCAACT-1 | 6616 | 3058 | 2 | PSE | 6.6928352 | 5.7139691 |
| 3 | APL_AAACCCACAGTTGAAA-1 | APL | AAACCCACAGTTGAAA-1 | 12781 | 3669 | 1 | MSE | 0.9346761 | 0.6447128 |
| 4 | APL_AAACCCAGTGCCTTTC-1 | APL | AAACCCAGTGCCTTTC-1 | 7138 | 3100 | 4 | PPP | -5.9457399 | 1.5383996 |
| 5 | APL_AAACCCAGTTTCGACA-1 | APL | AAACCCAGTTTCGACA-1 | 7364 | 2529 | 6 | PSE | 5.8819741 | 8.5710071 |
| 6 | APL_AAACCCATCTGCGAGC-1 | APL | AAACCCATCTGCGAGC-1 | 7114 | 2869 | 3 | CC | -0.7623315 | -0.4224323 |
rownames(meta) <- meta$cell_id
nrow(meta)
sce
class: SingleCellExperiment dim: 17396 10204 metadata(4): nucgenes genevar hvgs hvgs_vst assays(4): counts logcounts vst logvst rownames(17396): AT1G01010 AT1G01020 ... AT5G67630 AT5G67640 rowData names(10): ID chrom ... go_mitosis go_phloem_dev colnames(10204): APL_AAACCCAAGTAACGTA-1 APL_AAACCCACAACCAACT-1 ... sAPL_TTCGGTCAGCGACAGT-1 sAPL_TTTACCATCCTCTGCA-1 colData names(38): Sample Barcode ... slingPseudotime_4 slingPseudotime_5 reducedDimNames(36): MNN_logcounts MNN_logvst ... DIFFMAP_MNN_logcounts DIFFMAP_MNN_logvst mainExpName: NULL altExpNames(0):
seu <- as.Seurat(sce)
Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from MNN_logcounts_ to MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from MNN_logvst_ to MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from PC to PC_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to PC_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from PC to PC_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to PC_” Warning message: “Cannot add objects with duplicate keys (offending key: PC_), setting key to 'pca_logvst_'” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_logcounts_ to UMAP7logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_logvst_ to UMAP7logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_MNN_logcounts_ to UMAP7MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_MNN_logvst_ to UMAP7MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_logcounts_ to UMAP15logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_logvst_ to UMAP15logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_MNN_logcounts_ to UMAP15MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_MNN_logvst_ to UMAP15MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_logcounts_ to UMAP30logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_logvst_ to UMAP30logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_MNN_logcounts_ to UMAP30MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_MNN_logvst_ to UMAP30MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_logcounts_ to UMAP100logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_logvst_ to UMAP100logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_MNN_logcounts_ to UMAP100MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_MNN_logvst_ to UMAP100MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_logcounts_ to TSNE15logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_logvst_ to TSNE15logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_MNN_logcounts_ to TSNE15MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_MNN_logvst_ to TSNE15MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_logcounts_ to TSNE30logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_logvst_ to TSNE30logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_MNN_logcounts_ to TSNE30MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_MNN_logvst_ to TSNE30MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_logcounts_ to TSNE60logcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60logcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_logvst_ to TSNE60logvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60logvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_MNN_logcounts_ to TSNE60MNNlogcounts_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60MNNlogcounts_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_MNN_logvst_ to TSNE60MNNlogvst_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60MNNlogvst_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_” Warning message: “Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_logvst_'” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_” Warning message: “Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_mnn_logcounts_'” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_” Warning message: “All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_” Warning message: “Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_mnn_logvst_'”
seu
An object of class Seurat 17396 features across 10204 samples within 1 assay Active assay: originalexp (17396 features, 0 variable features) 36 dimensional reductions calculated: MNN_logcounts, MNN_logvst, PCA_logcounts, PCA_logvst, UMAP7_logcounts, UMAP7_logvst, UMAP7_MNN_logcounts, UMAP7_MNN_logvst, UMAP15_logcounts, UMAP15_logvst, UMAP15_MNN_logcounts, UMAP15_MNN_logvst, UMAP30_logcounts, UMAP30_logvst, UMAP30_MNN_logcounts, UMAP30_MNN_logvst, UMAP100_logcounts, UMAP100_logvst, UMAP100_MNN_logcounts, UMAP100_MNN_logvst, TSNE15_logcounts, TSNE15_logvst, TSNE15_MNN_logcounts, TSNE15_MNN_logvst, TSNE30_logcounts, TSNE30_logvst, TSNE30_MNN_logcounts, TSNE30_MNN_logvst, TSNE60_logcounts, TSNE60_logvst, TSNE60_MNN_logcounts, TSNE60_MNN_logvst, DIFFMAP_logcounts, DIFFMAP_logvst, DIFFMAP_MNN_logcounts, DIFFMAP_MNN_logvst
saveRDS(seu,"./tricycle/phloem_atlas_seu4.rds")
seu <- CreateSeuratObject(seu@assays$originalexp@counts, meta.data=meta)
seu <- NormalizeData(seu)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = nrow(seu))
all.genes <- rownames(seu)
seu <- ScaleData(seu, features = all.genes)
seu <- RunPCA(seu,features = VariableFeatures(object = seu))
Centering and scaling data matrix PC_ 1 Positive: AT5G16130, AT1G27400, AT2G37190, AT3G45030, AT5G27850, AT4G29410, AT2G34480, AT3G16080, AT3G62870, AT5G22440 AT1G22780, AT2G44120, AT3G11940, AT5G15200, AT1G09690, AT3G61110, AT5G20290, AT4G36130, AT2G39460, AT2G01250 AT3G22230, AT5G59850, AT3G07110, AT5G60670, AT5G48760, AT1G69620, AT5G52650, AT5G39740, AT1G09590, AT5G45775 Negative: AT1G76180, AT1G62480, AT5G65020, AT1G65930, AT2G02130, AT1G49240, AT5G56540, AT2G19760, AT3G21770, AT3G13520 AT1G10682, AT2G45820, AT1G20440, AT3G14310, AT1G05850, AT3G53420, AT5G17920, AT3G05890, AT1G12080, AT2G44790 AT2G30870, AT5G11970, AT3G23820, AT3G11930, AT3G56240, AT5G28050, AT1G35720, AT2G47730, AT1G55330, AT4G26320 PC_ 2 Positive: AT1G29470, AT3G58120, AT5G23860, AT2G45070, AT1G23490, AT5G14030, AT2G27030, AT1G70490, AT2G47170, AT5G15350 AT3G12610, AT5G48580, AT1G78040, AT1G04680, AT1G04820, AT3G03160, AT1G10630, AT2G27130, AT2G21160, AT3G43810 AT1G50010, AT5G17190, AT4G14960, AT4G24920, AT3G54810, AT2G19460, AT4G32460, AT5G03345, AT4G24190, AT3G25220 Negative: AT4G14010, AT2G37040, AT2G40890, AT4G34600, AT3G21240, AT5G57685, AT1G32100, AT5G64260, AT2G44300, AT5G54160 AT4G34050, AT5G48930, AT3G27890, AT1G44800, AT2G36830, AT2G30860, AT1G27030, AT2G19110, AT2G40900, AT3G23090 AT1G73660, AT1G18140, AT2G30490, AT3G26520, AT1G48320, AT3G10340, AT1G16350, AT1G06550, AT4G15340, AT1G48750 PC_ 3 Positive: AT3G28860, AT5G19780, AT4G19410, AT1G15690, AT1G09640, AT5G60390, AT5G02500, AT3G19820, AT1G07940, AT2G31360 AT4G13940, AT2G36530, AT4G38680, AT2G14890, AT4G34200, AT5G08690, AT5G14040, AT5G44340, AT3G20050, AT3G11830 AT5G62890, AT5G08670, AT3G04120, AT2G34250, AT1G04430, AT1G20330, AT3G13460, AT1G24510, AT5G20890, AT5G46290 Negative: AT4G24960, AT1G62045, AT3G17730, AT5G18130, AT2G35585, AT5G14110, AT1G69970, AT3G49360, AT5G50290, AT3G54040 AT5G02140, AT3G01940, AT1G33475, AT2G18380, AT1G18210, AT2G27740, AT1G06490, AT4G17510, AT1G65910, AT3G02677 AT3G18670, AT1G63310, AT1G11570, AT1G75720, AT3G08690, AT4G05220, AT1G05610, AT1G31935, AT3G03170, AT5G01070 PC_ 4 Positive: AT1G22840, AT4G32470, AT1G31812, AT1G72020, AT1G51650, AT4G21105, AT2G31490, AT5G64350, AT4G29480, AT1G15120 AT3G46430, AT5G47030, AT5G53560, AT4G20150, AT1G27450, AT3G52730, AT4G16450, AT5G08060, AT3G52300, AT4G34700 AT5G47890, AT4G37830, AT1G05205, AT5G59613, AT4G30010, AT1G76200, AT1G14450, AT3G10860, AT4G09030, AT5G05370 Negative: AT3G55760, AT3G24240, AT1G10760, AT5G14640, AT4G30020, AT2G37590, AT1G31880, AT1G06490, AT5G19090, AT1G65910 AT1G55760, AT1G54030, AT2G17840, AT5G61960, AT3G02260, AT4G05430, AT3G54040, AT5G01890, AT3G49360, AT5G60200 AT1G63300, AT2G18380, AT1G64390, AT4G39400, AT1G05200, AT1G05610, AT3G41768, AT1G72700, AT5G14110, AT2G37080 PC_ 5 Positive: AT5G26260, AT3G16420, AT1G54000, AT2G41800, AT2G47140, AT3G53480, AT2G43610, AT3G13790, AT5G60950, AT5G10130 AT3G15950, AT3G09260, AT5G02260, AT3G20370, AT1G53180, AT3G16460, AT2G04160, AT5G26280, AT2G17280, AT2G43535 AT2G30340, AT1G66270, AT1G78340, AT5G35940, AT1G54010, AT5G01040, AT5G14730, AT3G29034, AT4G23590, AT2G45400 Negative: AT2G26730, AT4G20270, AT1G68430, AT1G60390, AT3G15680, AT5G60200, AT2G16850, AT4G22120, AT1G72150, AT2G45470 AT3G52490, AT1G80690, AT5G25490, AT3G04670, AT5G57110, AT1G73850, AT3G63300, AT1G05470, AT2G18380, AT1G65910 AT3G28455, AT2G01910, AT5G54660, AT4G36620, AT1G70670, AT1G56230, AT1G63300, AT2G44830, AT3G45610, AT3G59680
seu <- FindNeighbors(seu, dims = 1:10)
seu <- FindClusters(seu, resolution = 0.5)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 10204 Number of edges: 303357 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9067 Number of communities: 14 Elapsed time: 1 seconds
seu <- RunUMAP(seu, dims = 1:10)
10:59:20 UMAP embedding parameters a = 0.9922 b = 1.112 10:59:21 Read 10204 rows and found 10 numeric columns 10:59:21 Using Annoy for neighbor search, n_neighbors = 30 10:59:21 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 10:59:22 Writing NN index file to temp file /tmp/RtmpSzDebc/file16ca8b2cbb33de 10:59:22 Searching Annoy index using 1 thread, search_k = 3000 10:59:25 Annoy recall = 100% 10:59:27 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 10:59:30 Initializing from normalized Laplacian + noise (using RSpectra) 10:59:30 Commencing optimization for 200 epochs, with 388778 positive edges 10:59:37 Optimization finished
head(seu@meta.data)
| orig.ident | nCount_RNA | nFeature_RNA | cell_id | sample | barcode | total_counts | detected_genes | cluster | annotation | UMAP1 | UMAP2 | RNA_snn_res.0.5 | seurat_clusters | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <fct> | <dbl> | <int> | <chr> | <chr> | <chr> | <int> | <int> | <int> | <chr> | <dbl> | <dbl> | <fct> | <fct> | |
| APL_AAACCCAAGTAACGTA-1 | APL | 35383 | 5176 | APL_AAACCCAAGTAACGTA-1 | APL | AAACCCAAGTAACGTA-1 | 35600 | 5210 | 9 | early/undetermined | 2.2210797 | -3.4197492 | 8 | 8 |
| APL_AAACCCACAACCAACT-1 | APL | 6456 | 3020 | APL_AAACCCACAACCAACT-1 | APL | AAACCCACAACCAACT-1 | 6616 | 3058 | 2 | PSE | 6.6928352 | 5.7139691 | 9 | 9 |
| APL_AAACCCACAGTTGAAA-1 | APL | 12731 | 3650 | APL_AAACCCACAGTTGAAA-1 | APL | AAACCCACAGTTGAAA-1 | 12781 | 3669 | 1 | MSE | 0.9346761 | 0.6447128 | 8 | 8 |
| APL_AAACCCAGTGCCTTTC-1 | APL | 6954 | 3072 | APL_AAACCCAGTGCCTTTC-1 | APL | AAACCCAGTGCCTTTC-1 | 7138 | 3100 | 4 | PPP | -5.9457399 | 1.5383996 | 4 | 4 |
| APL_AAACCCAGTTTCGACA-1 | APL | 7276 | 2483 | APL_AAACCCAGTTTCGACA-1 | APL | AAACCCAGTTTCGACA-1 | 7364 | 2529 | 6 | PSE | 5.8819741 | 8.5710071 | 11 | 11 |
| APL_AAACCCATCTGCGAGC-1 | APL | 7045 | 2849 | APL_AAACCCATCTGCGAGC-1 | APL | AAACCCATCTGCGAGC-1 | 7114 | 2869 | 3 | CC | -0.7623315 | -0.4224323 | 1 | 1 |
DimPlot(seu, reduction = "umap", group.by = "annotation")
sce <- as.SingleCellExperiment(seu)
## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)
The number of projection genes found in the new data is 242.
sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi
## Plot
print(plot_ccposition_den(sce$tricyclePosition,
sce$orig.ident, 'sample',
bw = 10, fig.title = "Kernel density of \u03b8") +
theme_bw(base_size = 14))
## Assign position to stage
tricycleCCStage <- sce$tricyclePosition/pi
G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
tricycleCCStage[G1_idx]="G1/G0"
tricycleCCStage[S_idx]="S"
tricycleCCStage[G2_idx]="G2M"
sce$tricycleCCStage <- tricycleCCStage
## Schwabe Stage
#sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
## Table
#seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
#print(table(seu$time.anno.Li, sce$tricycleCCStage))
## Prepare seu
#seu <- prep(seu)
#plot_anno(seu)
## Sce to seu
seu$tricyclePosition <- sce$tricyclePosition/pi
seu$tricycleCCStage <- sce$tricycleCCStage
#seu$SchwabeCCStage <- sce$CCStage
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
point.size = 4, point.alpha = 1) +
theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
DimPlot(seu, reduction = "umap", group.by = "tricycleCCStage")
saveRDS(seu,"./tricycle/phloem_atlas_seu4_simplified.rds")
# Set the identity of interest to cell type
Idents(rc.integrated) <- "tricycleCCStage"
# Use SCT assay
DefaultAssay(rc.integrated) <- "integrated"
de.list <- FindAllMarkers(rc.integrated,
max.cells.per.ident = 10000,
only.pos=T,
test.use="wilcox")
Calculating cluster G2M
de.list <- de.list %>% filter(p_val_adj < 0.01 & avg_log2FC > 0.25)
table(de.list$cluster)
G2M G1/G0 S 133 9 103
de.list %>% filter(cluster=="G1/G0")
| p_val | avg_log2FC | pct.1 | pct.2 | p_val_adj | cluster | gene | |
|---|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <fct> | <chr> | |
| AT5G02820 | 1.728994e-120 | 2.7087332 | 0.138 | 0.378 | 4.236034e-118 | G1/G0 | AT5G02820 |
| AT1G04710 | 2.012929e-90 | 2.8375044 | 0.120 | 0.281 | 4.931675e-88 | G1/G0 | AT1G04710 |
| AT2G03670 | 4.696481e-50 | 4.8883449 | 0.138 | 0.325 | 1.150638e-47 | G1/G0 | AT2G03670 |
| AT3G13060 | 2.300513e-25 | 1.9837473 | 0.133 | 0.147 | 5.636256e-23 | G1/G0 | AT3G13060 |
| AT2G42190 | 1.941702e-21 | 0.8564682 | 0.126 | 0.202 | 4.757170e-19 | G1/G0 | AT2G42190 |
| AT4G14980 | 3.124803e-08 | 2.9606000 | 0.111 | 0.129 | 7.655767e-06 | G1/G0 | AT4G14980 |
| AT1G05520 | 1.283976e-07 | 1.8709291 | 0.183 | 0.239 | 3.145741e-05 | G1/G0 | AT1G05520 |
| AT1G128001 | 1.191791e-06 | 0.7292746 | 0.141 | 0.264 | 2.919888e-04 | G1/G0 | AT1G12800 |
| AT5G63120 | 2.330971e-06 | 4.8604293 | 0.110 | 0.179 | 5.710880e-04 | G1/G0 | AT5G63120 |